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We present density functional perturbation theory for electric field perturbations and ultra-soft 
pseudopotentials. Applications to benzene and anthracene molecules and surfaces are reported as 
examples. We point out several issues concerning the evaluation of the polarizability of molecules 
and slabs with periodic boundary conditions. 
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I. INTRODUCTION 

Ultrasoft pseudo-potentials (US-PPs)^ are employed in 
large scale electronic structure calculations because they 
allow precise and efficient simulations of localized ?>d and 
2p electrons with plane- waves basis sets. Several ab-initio 
techniques have been adapted to US-PPs. Important 
examples comprise the Car-Parrinello molecular dynam- 
ics^ the Berry phase approach to the macroscopic polar- 
ization of solidS(^ the ballistic conductance of open quan- 
tum systems, i and also density functional perturbation 
theory (DFPT),^*S^ DFPTi2*a addresses the response of 
an inhomogeneous electron gas to external perturbations, 
giving access to the lattice dynamics, to the dielectric, 
and to the elastic properties of materialsfiS. DFPT for 
lattice dynamics with US-PPs has been presented in de- 
tail in Refs.00, whereas the treatment of an electric field 
perturbation and of the dielectric response has been only 
briefly sketched in the literaturei^ The purpose of this pa- 
per is to describe our implementation of DFPT for elec- 
tric field perturbations, to validate it with US-PPs and to 
apply it to some examples. We calculate the polarizabil- 
ities of two molecules, benzene and anthracene, and the 
dielectric properties of the (010) surface of benzene and of 
the (100) surface of anthracene. We validate our DFPT 
implementation with US-PPs in two ways. First we focus 
on the electronic density induced, at linear order, by an 
electric field. This density is calculated by DFPT and 
by a self-consistent finite electric field (FEF) approach. 
Second, we compare the polarizability inferred from the 
dipole moment of the induced charge and from the dielec- 
tric constant of the periodic solid simulated within the 
super-cell approach. In both cases we find a very good 
agreement between DFPT and FEF. 

Using plane- waves and pseudo-potentials, as imple- 
mented in present codes, it is not possible to study truly 
isolated molecules or surfaces. A super-cell geometry 
must be adopted by creating a fictitious lattice made of 
periodically repeated molecules or slabs (with two sur- 
faces) separated by enough vacuum and the convergence 
of the calculated properties must be checked against the 
increase of the size of the vacuum region. In this re- 
spect, to converge the dielectric properties is particularly 
difficult because, due to the long range of electrostatic 
interactions, periodic boundary conditions may generate 



a spurious electric field which become negligible only at 
very large vacuum spacing. In this paper we show that, 
accounting for this spurious electric field, it is possible to 
reduce considerably the vacuum needed to converge the 
polarizability both by DFPT and by FEF. 

We note in passing that in solid benzene and an- 
thracene there are two distinct chemical bonds: strong 
covalent bonds within a single molecule and weak bonds 
between molecules responsible for the cohesion of the 
solid. Within density functional theory, in the local den- 
sity approximation (LDA) or in the generalized gradi- 
ent approximation (GGA), one cannot account for the 
van der Waals forces which are important for molecule- 
molecule interactions and therefore the geometry of the 
system cannot be determined from energy minimiza- 
iioum^i^ However, once the geometry has been adjusted, 
the other calculated properties often turns out to be rea- 
sonably described by LDA or GGA. In this paper, while 
we allow the geometry of the isolated molecules to be fully 
relaxed in order to minimize energy, we borrow from ex- 
perimentii the orientations of the molecules and the cell 
sizes for bulk benzene and anthracene. The geometry of 
the surfaces is assumed to be identical to the truncated 
bulk. 

The outline of the paper is the following: Section II 
contains an expression for the dielectric constant with 
US-PPs. Section III describes our implementation of the 
FEF approach and other technical details. Section IV 
is devoted to the study of the polarizability of the ben- 
zene and anthracene molecules, while in Section V the 
dielectric properties of slabs of benzene and anthracene 
are discussed. 



II. DIELECTRIC CONSTANT WITH 
ULTRASOFT PSEUDO-POTENTIALS 

The macroscopic dielectric tensor of an insulating solid 
is defined a.S'^ 



eafB = Safj + 47r 



(1) 



a and /? are Cartesian coordinates and is the deriva- 
tive of the electronic polarization with respect to the 
macroscopic screened electric field E. We neglect atomic 
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relaxations and therefore all our considerations regard 
the so called "clamped-ions" dielectric constant e^o where 
only the dielectric contribution of electrons is accounted 
for. 

In order to calculate , we begin with the derivative 
of the electronic density with respect to an electric field. 
With US-PPs, the electronic density is written as: 



(V,|^(r)|^,) 



(2) 



where the integral is over the volume of the solid, made 
up of iV unit cells of volume fi. Born- von Karman bound- 
ary conditions are assumed for the wave- functions. The 
electron charge is (— e). 

For convenience, we define the functions \4>") = 
J d^r eraK{r)\ipi), and introduce the projector into the 
empty states manifold Pc = J2c IV'c)('0c|<S', where S is the 
overlap matrix (see below). With these definitions, Eq.[Sl 
becomes: 



where the sum runs over the occupied states and the 
kernel /ir(r;ri,r2) is: 

i4:(r;ri,r2) = S{r - ri)S{r - + 
^ Qlg\r - RiM^'Hr, - Ili)f3*J'Hv, Rj). (3) 
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The augmentation functions QZm\'<^) and the projector 
functions f3n'"^\r) are calculated togheter with the PP 
and are localized about the atoms at position R/ji 

The electronic charge linearly induced by an electric 
field is therefore: 

c.c. means complex conjugate. 
Using Eq. we get as: 



dEa 



2e 



Y.jd\{^\r^K{v)\^,)+c.c., (5) 



dEp Nn^^dEp^^^^'^^ 



(6) 



The functions \(j)f) and hence are well defined in 

a finite system. In a periodic solid can be de- 

fined as done with norm conserving pseudo-potentials. 
We exploit the relation between the off-diagonal matrix 
elements of the r operator between non-degenerate Bloch 
states, and the matrix elements of the velocity operator. 
With US-PPs, we have: 



{■)jjj\Sra\->Jji) 



{ipj \ [H - eiS,ra] iV'i) 



£'i - e-i 



(7) 



where H is the unperturbed Hamiltonian and £i are the 
unperturbed eigenvalues. The overlap matrix S isi 



5(ri,r2) = 5(ri - r2) + ^ Q^iil^ f^^J'Hr, - Ki)f3*^^'Hr2 ~ R/), (8) 

Inm 

where the coefficients qZm^ — J d^r QZini'r) are the integrals of the augmentation functions. 
Using Eq. 13 we see that Pc^'alV'i) are the solutions of the linear system: 

[H + Q- e^S] PcraH^) = Pi [H - e^S, ra] \^^), (9) 

where both the left and the right hand sides are lattice periodic. Q is added to the linear system in order to make 
it non singular as explained in detail in Ref. (See Eq. 30 and Eq. 72). With norm conserving pseudopotentials, 
in insulators, Q is proportional to the valence band projector. Its generalization to US-PPs is given in Ref. ifl (see 
discussion after Eq. 29). By solving this Hnear system, we obtain Pcra\fpi), while we need P^ralfpi) to compute 
Pj|0f ). Since SPc = PjS", we have S'Pcr„|V'») PjS'ralV',), hence the functions Pj|0f ) are: 

plm = sPcerM - Pi E + Pi E ifuMiP'M^). (10) 

Inva Inm 

where we defined the integral /f„„ = / d^r eraQZini^ — R/). 

To proceed further, we need the first order derivative of the electronic wave- functions {ipi) with respect to an electric 
field. We can calculate these quantities to linear order in perturbation theory. The overlap matrix S does not depend 
on the electric field while the differentiation of the Kohn and Sham potential yields: 



dVKS 
dEp 



d-^r 



, dVHxc{r) 



K{r), (11) 
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where Vhxc is the Hartree and exchange and correlation potential. Applying to Eq. (19) of Ref. 0, we obtain: 



[H + Q~ e,S] Pc 



dipt 



= -P} 



dEa 



X(r)|V. 



(12) 



The self-consistent solutions of this linear system, to- 
gether with are substituted into Eq.^lto calculate 
4^. Finally, the dielectric tensor is calculated via Eq.^ 



III. TECHNICAL DETAILS 

We validate the theory by comparison with a self- 
consistent FEE method. A FEE is simulated as sug- 
gested by Kunc and Restat^^ii^ a sawtooth-like poten- 
tial is added to the bare ionic potential. This method, 
when applied to periodic systems, is not as efficient as 
other recent approachesiSiiLiSii^iSS because it requires 
the simulation of a super-cell, but its implementation is 
simple. Our systems, molecules and slabs, require al- 
ready a super-cell and the method of Kunc and Resta 
suits adequately our goals. The molecules and the slabs 
are inserted in the region where the sawtooth-like poten- 
tial is linear like the potential of an electric field (E^z): 
4>(r) = — Esi -r. To ensure periodicity, the slope of the po- 
tential is reversed in a small region in the middle of the 
vacuum. For finite super-cells and small enough fields, 
the electrons are only slightly polarized by the field, the 
systems maintain a well defined energy gap between oc- 
cupied and empty electronic states and electrons do not 
escape into the vacuum. In these simulations, the macro- 
scopic electric field E is zero because, solving the Pois- 
son equation we assume, as a boundary condition, the 
periodicity of the Hartree potential. The sawtooth-like 
potential describes a microscopic electric field Ejj which 
vanishes averaging over macroscopic distances. On our 
finite systems, this microscopic field has the same effect 
of a macroscopic electric field. Indeed, in DEFT, as de- 
rived in previous Section, the perturbation is a genuine 
macroscopic electric field E and, as shown in the follow- 
ing, the electronic density induced, at linear order, by 
the field Esz and the electronic density induced by the 
macroscopic electric field E are equal to each other when 
Esi = E. 

Both the FEE approach and DEPT for US-PPs have 
been implemented in the PWscf'^^ package. All calcu- 
lations are carried out within the GGA approximation. 
The expression of the exchange and correlation energy 
introduced by Perdew, Burke and Ernzerhof (PBE)S 
is used in the GGA functional. The US-PPs of hydro- 
gen and carbon have the parameters described in Ref. 0. 
Plane waves up to an energy cutoff of 30 Ry for the wave- 
functions and 180 Ry for the electron density, are used. 
Only the F point is used for the Brillouin zone (BZ) sam- 
pling in the molecular calculation while a 3 x 3 Monkhorst 
and Pack mesh^'^ of k-points is used for sampling the two- 
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FIG. 1: Benzene and anthracene molecules. Distances and 
angles are reported in Tables Q] and |n] for benzene and an- 
thracene respectively. The polarizability tensors are given 
with respect to the axes shown in the figure. 



dimensional BZ of the slabs. 



IV. MOLECULES 

Benzene (CgHg) and anthracene (C14II10) are planar 
molecules, their geometries are shown in Fig. ^ We op- 
timize the geometries and our theoretical bond lengths 
and angles are reported in the Tabs. HI and |nl and 
compared with previous calculations^^ and with exper- 
iment.^^ Overall, the agreement between theory and ex- 
periment is good and our values compare well also with 
the more precise B3LYP resultsi24 

As a first test of DEPT, we calculate the electronic 
density A?! ~ induced at linear order by an electric 
field. Fig. 12 (Fig. EJ shows the density induced in ben- 

this work Ref. 24 Ref. 1^ 
(PBE) (B3LYP) (expt.) 

C-C 1.396 1.399 1.398 (n), 1.392 (x) 

C -H 1.091 1.092 1.090 (n) 

TABLE I: Theoretical GGA optimized geometry of the ben- 
zene molecule compared with previous theoretical worh— 
(B3LYP with localized basis set) and experiment^ Bond 
lengths are in A. The symbols are defined in Fig. Ab- 
breviations: (n) neutron, (x) x-ray diffraction. 
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FIG. 2: Electronic density induced by an electric field either parallel (a) or perpendicular (b) to the plane of benzene calculated 
by DFPT. (1) Plane perpendicular to the z axis passing through (0.0,0.0, —2.0) a.u.; (2) plane xz. Contours are in units of 
10~^ elec./(a.u.)^. The contours correspond to a field 0.5 a.u. (2.57 x 10^ V/cm). 



zene (anthracene) by an electric field either parallel or 
perpendicular to the molecular plane. The same An has 
been calculated by FEF via a numerical differentiation 
of the self-consistent density with lE^il = 10"'^ a.u. and 
|Es;| = a.u.. In all cases, to the scale of the figures, 
the An calculated by DFPT and the An calculated by 
FEF coincide. In Fig. 0] we report, as an example, the 
difference between two An to an enlarged scale. This 
error is due to nonlinear effects (present in the FEF re- 
sults but not in DFPT) as well as to numerical noise. Its 
magnitude, lower than 1% of An, is similar with norm- 
conserving PPs. 

Now we can address the molecular polarizabilities. The 
polarizability of a molecule is a tensor which measures, 
at linear order, the tendency of the molecule to change 
its dipole moment when inserted into an electric field. It 
is defined by the relationship — J2p '^a/3'Eiioc,0, be- 
tween the dipole moment p of the induced charge den- 
sity, and Eioc the electric field acting on the molecule. In 
general, the polarizability tensor a can be made diago- 
nal in the principal axes. The point group of benzene. 



Dqii, and of anthracene, D2h, are centrosymmetric, hence 
these two molecules have no permanent dipole moment. 
Their principal axes are shown in Fig.^ In these axes, a 
has two independent components in benzene {a^x = ctyy 
parallel and azz perpendicular to the molecular plane), 
and three independent components in anthracene {axx, 
Oiyy , Oizz , parallel to the short molecular axis, parallel to 
the long molecular axis and perpendicular to the molec- 
ular plane, respectively). We first extract the polariz- 
abilities from the dipole moment of the induced charge 
density — eAn. As shown in Figs. [21 and |21 the An in- 
duced in our molecules are localized and well separated 
by their periodic images, so that p can be calculated 
by a numerical integration (p = — e / An(r)r cPr) over 
the super-cell volume. The small differences between in- 
duced charges calculated by DFPT and by FEF lead to 
differences smaller than 1% in the dipole moments (for 
instance, in benzene, for L ~ 2A a.u. and E = 0.5 a.u., 
Px = 42.79 a.u. with DFPT and px = 42.85 a.u. with 
FEF). The polarizabilities differ also by less than 1% be- 
cause, at linear order, E;oc acting on the molecules is 
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FIG. 3: Electronic density induced by an electric field either parallel (a,b) or perpendicular (c) to the plane of anthracene 
calculated by DFPT. In (a) the field is parallel to the short axis, in (b) it is parallel to the long axis of the molecule. The plane 
is perpendicular to the z axis passing through the point (0.0,0.0, —2.0) a.u.. Contours are in units of 10"^ elec./(a.u.)^. The 
contours correspond to a field 0.5 a.u.. 



the same. Indeed, the local field can be estimated rec- 
ognizing that the periodic solid simulated in the super- 
cell approach, has a macroscopic polarization P = p/fi. 
Hence, the local field which acts on the molecules is 
given by the Lorentz formula E/oc = Es; -t- (FEF) 
or E,oc = E + ^P (DFPT). Of course, these formulas 
are valid for isotropic solids and cubic super-cells, how- 



bond 


this work 


Ref. J4 




Ref.J^ 


(PBE) 


(B3LYP) 


(expt.) 


(expt. solid) 


Ci -Cs 


1.398 


1.403 


1.392 


1.403 


Ca "Ct 


1.425 


1.432 


1.437 


1.445 


C^ — Cll 


1.370 


1.372 


1.397 


1.374 


Cs — C4 


1.446 


1.447 


1.437 


1.425 


Cll — C12 


1.421 


1.428 


1.422 


1.412 


Ci^m 


1.092 


1.094 




1.121 


Cll — H7 


1.090 


1.092 




1.139 


Cj — H-j. 


1.091 


1.093 




1.153 


angle 


C3C1C5 


121.74 


121.80 




120.37 


C3C7C11 


120.89 


120.97 




120.31 


H3C7C11 


120.55 


120.53 




123.48 


H7C11C12 


119.34 


119.43 




124.93 



TABLE II: Theoretical GGA optimized geometry of the an- 
thracene molecule compared with previous theoretical work— 
(B3LYP with localized basis set) and experiment!— The ex- 
perimental geometry of solid anthracene is also reported for 
comparison.^ Bond lengths are in A, angles in degrees. The 
symbols are defined in Fig. Q 



ever the isotropic formula is sufficient to correct the main 
effects of the periodic boundary conditions when the elec- 
tric field is parallel to one principal axis and the super- 
cell is cubic. We evaluate the dipole moment of benzene 
and anthracene for several sizes of the box. In Tab. IIIII 
we report the polarizabilities calculated with or without 
the Lorentz corrections. In benzene, with a super-cell of 
20 a.u. and a polarizability axx = 83.6 a.u., the local 
field is 4.4% larger than Egi or E. In our largest super- 
cell {L = 32 a.u.) the local field is still 1% larger than 
Esi or E. From Tab. IIIII we can see that the inclusion 
of the Lorentz correction increases the convergence rate 
of axx and azz- A box of 20 a.u. is sufficient to give 
values converged within 1%. The anthracene molecule 
is about 18 a.u. long in the y direction and therefore 
the molecules can be considered as truly isolated only for 
super-cell sizes larger than 28 a.u.. At smaller box sizes, 
not only electrostatic, but also direct molecule-molecule 
interactions are important. At L = 28 a.u., including the 
Lorentz correction, axx and azz are already converged 
within 1%. ayy is instead more difficult to converge. A 
super-cell of about 50 a.u. is needed to reduce the local 
field effects below 1%. Instead, including Lorentz correc- 
tions, a box size of 32 a.u. is sufficient to converge ayy 
within 1%. 



Besides the direct approach, molecular polarizabilities 
can be evaluated also starting from the dielectric con- 
stant (Eq. of the periodic solid simulated in the super- 
cell approach. The anisotropic Clausius-Mossotti for- 
mula which derives from the approximate Lorentz field 
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Benzene Anthracene 



L (a.u.) 




Olzz 








20 


83.6 (87.4) 


44.7 (45.7) 








— 


24 


83.5 (85.7) 


44.9 (45.5) 


170 (179) 


317 (351) 


85.5 (87.8) 


28 


83.5 (84.8) 


44.9 (45.3) 


171 (177) 


313 (333) 


86.3 (87.8) 


30 


83.5 (84.6) 


45.0 (45.3) 


172 (176) 


310 (327) 


86.6 (87.9) 


32 


83.4 (84.3) 


45.0 (45.2) 


172 (176) 


309 (322) 


86.8 (87.8) 


50 






172 (173) 


306 (309) 


86.5 (86.8) 



TABLE III: Independent components of the molecular polarizability tensor of benzene and anthracene (in a.u.) calculated from 
the dipole moment of — eAn including Lorentz field corrections. Values in parenthesis are calculated neglecting Lorentz field 
corrections. 

Benzene Anthracene 

L (a.u.) e^x {Clxx) £zz (Qzz) <^xx [a^x) ^yy [oi-yy) ^zz (Qzz) 



20 1.1372 (83.5) 1.0716 (44.5) — — — 

24 1.0778 (83.4) 1.0412 (44.7) 1.1625 (170) 1.3329 (330) 1.0793 (85.0) 

28 1.0485 (83.3) 1.0259 (44.8) 1.1013 (171) 1.1908 (313) 1.0501 (86.1) 

30 1.0393 (83.3) 1.0210 (44.8) 1.0819 (171) 1.1519 (311) 1.0407 (86.3) 

32 1.0323 (83.3) 1.0173 (44.9) 1.0673 (171) 1.1233 (309) 1.0335 (86.5) 



TABLE IV: Dielectric constant of benzene and anthracene molecules in a cubic super-cell as a function of the box size. In 
parenthesis we report the polarizability (in a.u.) obtained from Eg. 1131 



given above is 



,26 



SURFACES 



47r eaa + 2 



(13) 



In Tab. II VI we report the dielectric constant as a function 
of the super-cell size and the polarizabilities calculated 
via Eq. El for both molecules. The convergence rate 
of the polarizabihtics is similar to the convergence rate 
found in Tab. IIIII including the Lorentz field. In all cases 
the difference of the final polarizabilities in Tab. IIIII and 
in Tab.Cvlis below 1%. 

The final calculated components of the benzene polar- 
izabilities are axx — 83.4 a.u. and azz = 45.0 a.u.. These 
values give a mean polarizability a = ^{2axx + cizz) = 
70.6 a.u. and are in good agreement with the results 
of previous theoretical works and with experiment. We 
report in Tab. previous theoretical data and experi- 
mental values. For anthracene we get axx = 171 a.u., 
ayy = 306 a.u., azz — 86.8 a.u.. These values give a 
mean polarizability a = 188 a.u.. Experimental values 
of the anisotropic components of the polarizability of an- 
thracene exist for both diluted benzene solutions and for 
solid anthracene. The values of axx ranges from 139 a.u. 
to 174 between 238 a.u. and 292 a.u. and azz 

between 80 a.u. and 115 a.u. (see Tab. 0). Note that 
the ionic contribution to the polarizability, which is not 
accounted for in our calculation, is present in the experi- 
ment of Ref . but not in those of Refs. 36 3l However 
it has been estimated that the ionic contribution is only 
of the order of 5%iSS 



Some surfaces of molecular crystals can be cleaved by 
cutting only weak molecule-molecule bonds. Such sur- 
faces are expected to remain insulating and to have di- 
electric properties similar to the bulk. In this section, we 
present the dielectric properties of two of these surfaces: 
the (010) surface of benzene and the (100) surface of an- 
thracene. Moreover, we show how to calculate the dielec- 
tric properties of a single slab from the dielectric constant 
of the solid simulated in the super-cell approach. 

The (010) benzene surface is simulated by a six-layers 
slab and two molecules per layer (see Fig.lSJt)- The super- 
cell is orthorhombic with sizes 13.78 a.u. xL a.u. x 12.74 
a.u., where L depends on the vacuum between slabs. The 
(100) surface of anthracene is described by a four-layers 
slab and one molecule per layer (see Fig.Eb). The super- 
cell is monoclinic with sizes |a| = L/siwy, |b| = 11.41 a.u. 
and |c| = 21.14 a.u., where the angle between a and c is 
7 = 124.7°. The long axis of the molecules is approxi- 
mately parallel to the c vector. Fig. shows the shape 
of the sawtooth-like potential used in the FEF simula- 
tions. The applied field is normal to the surface, in the 
direction (010) in benzene and in the direction of b x c 
in anthracene. The band structures of these slabs have 
been reported in Ref. [s^- 

We begin by a comparison of the induced electronic 
density calculated by DFPT and by FEF. In order to 
visualize the induced density, we introduce the planar 
average of An, defined as Ana.„(A) = 1/S An(r) dS. 
Here A is a coordinate along the surface normal and S 
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Q 


Method /reference 


Anthracene 






164 


287 


81.2 


177 


APSC + B3LYP/6-311G(d,p)^ 


158 


266 


81.8 


167 


APSC + MP2 corrections^ 


172 


303 


86.5 


187 


APSC + TZVP-FIp2i 


171 


306 


86.8 


188 


this work (GGA-PBE) 


191 


337 


96.2 


208 


DFT with TZVP-FIP basis^ 


169 


240 


105 


171 


EFNMR - BS (static)^ 


173 


238 


103 


171 


KE, EP, CME - BS (dynamic)^ 


165 


242 


107 


171 


KE, EP, EST - BS (dynamic)^ 


174 


272 


80.3 


175 


CR (dynamic)^ 


139 


292 


115 


182 


CR (dynamic)^ 


Benzene 








74.2 




39.5 


62.6 


HF^^ 


76.0 




41.3 


64.4 


HF + MP2Si 


76.8 




37.2 


63.6 


HF with DZP' basisffi 


78.2 




37.2 


64.5 


as above + MP'M 


79.5 




45.2 


68.1 


HF with POL basisffi 


81.6 




45.2 


69.5 


as above + MPOiS 


83.4 




45.0 


70.6 


this work (GGA-PBE) 


83.7 




44.9 


70.9 


LDA, Gaussian basis^^ 


85.0 




45.6 


72.2 


as above - different basis— 


74.9 




49.9 


66.6 


Ref. 37 


74.9 




50.6 


66.8 


KE, EP, CME (dynamic) ^ 



TABLE V: Theoretical and experimental values of benzene 
and anthracene polarizabilities. All values are in atomic units. 
Abbreviations: Atomic polarizability in a self-consistent lo- 
cal field (APSC), benzene solution (BS), Kerr effect (KE), 
Cotton- Mouton effect (CME), electron polarization (EP), 
crystal refraction (CR), empirical estimate (EST), Hartree- 
Fock (HF), M0ller-Plesset correction of 2nd order (MP2). 



Benzene Anthracene 



FIG. 4: Density profiles along the line (A,0.0,-2.0) a.u., -12 

a.u.< A <12 a.u., (dashed line) and the difference between L (a.u.) a (Eq.[T6|l a L (a.u.) a (Eq.[T6l a 

DFPT and FEE (continuous line). Upper panel for benzene 65.0 500.7 1116 36.17 298.6 523.9 

in a cubic super-cell with size 24 a.u. The lower panel for 53 g 500.9 1044 39.46 298.0 491.3 

anthracene in an orthorhombic super-cell with sizes 24 a.u. x r,g g ^qq g ggg g ^2 75 298 1 468 2 



36 a.u. X 15 a.u. 



80.0 500.9 907.9 46.03 298.1 449.9 
85.0 500.9 866.5 49.32 298.1 435.1 



is the area of one surface unit celL The integration is 
on cross sections parallel to the surface. In Fig. |H1 we 
report Auav calculated by DFPT and by FEF. In the 
latter case the numerical differentiation is done taking 
the field Egi = 10~^ a.u. and Egi — a.u. in benzene 
and Esi = 2 x 10^'^ a.u. and E^i — a.u. in anthracene. 
To the scale of the figure the planar averages calculated 
by FEF or by DFPT coincide. The error, reported to an 
enlarged scale in the insets, is always lower than 1% of 
the value of Auav 

Now we move on to the dielectric properties of the 
slabs. We begin by extracting the polarizability of the 
slab from the induced dipole moment per unit surface. 
We define the polarizability of a slab a, starting from 



TABLE VI: Polarizability of the benzene and anthracene slabs 
(in a.u.) as a function of the size of the super-cell. Values 
corrected for the electric field due to the periodic boundary 
conditions and uncorrected values are both reported. 

the relationship: 

m = ^Eioc, (14) 

where m = — e J^^AnQ^(A)A dX is the induced dipole 
moment per unit surface area and Eioc is the electric 
field, perpendicular to the surface, which induces the 
dipole moment on the slab {Eioc is the uniform electric 
field that remains in the vacuum between two slabs af- 
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(a) 




(b) 




//// A 



FIG. 5: Geometry of the benzene (a) and anthracene (b) 
slabs. The unit ceU is shown by a black frame. In (a) the 
applied sawtooth-like potential is also indicated. 



ter subtraction of all short range inhomogeneous fields). 
In a periodic slab geometry, the dipole moment of the 
slab is the origin of a sizable electric field. Actually, as 
described in Refs. ^^32, in an isolated slab a dipole mo- 
ment per unit surface area produces an the electrostatic 
potential which has different constant values in the vac- 
uum on the left and on the right part of the slab. In 
a repeated slab geometry, this jump cannot be accomo- 
dated with periodic boundary conditions. Therefore, an 
artificial uniform electric field appears in the super-cell 
in order to restore the periodicity of the electrostatic po- 
tential. Before applying Eq.^l this field has to be added 
to Esi or to E in order to calculate Eioc which ac tually 
induces the dipole moment on the slab. As in Refs. l3ll32 
we evaluate Eioc as: 



El, 



Eo 



IT' 



(15) 



where L is the length of the super-cell in the direction 
perpendicular to the surface and either Eq = Egi (FEF) 
or Eq= E (DFPT). Using Eqs. [Tland[TSl we calculate 
the polarizability of a slab from the dipole moment per 
unit surface area as: 



a 



4tt1 + x 



(16) 



3 
o 

0) 



.« -0.01 




Z2 
CD 



< 




0.0 



10.0 20.0 
X,(a.u.) 



30.0 



FIG. 6: Planar average Auav of the electron density induced 
by an electric field on the benzene slab (upper panel) and 
anthracene slab (lower panel). The difference between DFPT 
and FEF is shown in the insets. An corresponds to an electric 
field of 0.5 a.u.. Vertical lines indicate the position of the 
centers of the molecules in each layer. 



where x = ATrm/ (LEq). In Tab. IVll we report a calcu- 
lated by Eq. 1161 and compare with a obtained without 
including the field due to the periodic boundary condi- 
tions (a = Sm/Eo). It is found that the polarizability 
calculated by Eq.^|is already converged at the smallest 
vacuum space for both benzene and anthracene. Tab. I VII 
shows that the field due to the periodic boundary condi- 
tions has the same magnitude of the applied electric field. 
For instance, in benzene (anthracene), at the minimum 
slab-slab distance, for L = 65.0 a.u. (L = 36.2 a.u.), Eioc 
is about 123% (75%) larger than Eq. A vacuum size of 
about 3530 a.u. (1530 a.u.) would be necessary to re- 
duce the effect of the field due to the boundary conditions 
below 1%. 

Now, we can apply a similar argument and calculate 
the polarizability of a slab using the dielectric constant 
(Eq. of the periodic solid simulated in the super-cell 
approach. If n is a unit vector normal to the surface, the 
relevant dielectric constant for fields normal to the sur- 
face is e = X]a/3 ^apnanp. In benzene the surface normal 
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400 



200 




FIG. 7: Benzene slab polarizability as a function of the num- 
ber of layers. The data are interpolated with a linear fit (see 



is along the (010) direction and e — e 
it is in the xy plane and e = 



where n = 



yy, in anthracene 

_ ^^xy^x^y ~t~ ^yy^y 

0). We report in Tab.EUlthe dielec- 
tric constant as a function of the length of the super- 
cell. The dielectric constant depends on the size of the 
super-cell. Instead, the polarizability a becomes con- 
stant as soon as the vacuum size avoids the direct inter- 
action between neighbouring slabs. From Eq.^ we have 
e = l + ^aEioc/E and, by using Eq. ^1 we obtain the 
expression of the polarizability of the isolated slab as a 
function of the dielectric constant: 



a 



An 



1 



(17) 



We report in Tab. lVlJ the polarizability for each vacuum 
distance calculated by Eq. El These values converge to 
the same values of Tab. |^ and the convergence rate is 
similar. 

Now, we compare the dielectric properties of the sur- 
face with those of the bulk. Moreover, we extract the 
dielectric constant of the bulk (in the direction of the 
surface normal) from the slabs polarizabilities. We use a 
method inspired by Ref. .'i.'i. and restrict our attention to 
the (010) surface of benzene. For a fixed super-cell size 
{L = 68.9 a.u.) we calculate the polarizabihty of ben- 
zene slabs with different numbers of atomic layers. Fig.[7| 
shows a{N) calculated as described above for iV = 2 to 

= 6. It is found that a{N) increases linearly with the 
number of layers. We can understand this behavior using 
Eq. (5) of Ref. [s^- The polarization P/v of an isolated N 
layers slab is approximately 



N 



2Spo 
Nflf 



(18) 



where Coo is the surface charge of a very thick slab in 
which the surface contribution to the polarization is neg- 





Benzene 






Anthracene 




L (a.u.) 


e 


a (Eq.[T7l L (a.u.) 


e 


a (Ea.\m 


65.0 


2.2304 


500.9 


36.17 


1.7549 


298.7 


68.9 


2.0850 


500.9 


39.46 


1.6489 


298.1 


75.0 


1.9160 


500.9 


42.75 


1.5705 


298.1 


80.0 


1.8123 


500.9 


46.03 


1.5090 


298.1 


85.0 


1.7296 


500.9 


49.32 


1.4595 


298.1 



TABLE VII: Dielectric constant of the "solid" made up of 
periodically repeated benzene and anthracene slabs separated 
by vacuum space, as a function of the length (L) of the super- 
cell. The slab polarizability a (in a.u.) is evaluated using 
Eq.[T7| 



ligible. Poo is the sum of the surface-dipoles which ac- 
counts for the difference between the polarizability of 
the surface layers and that of the bulk, fis is the vol- 
ume of a bulk unit cell with two layers. In our exam- 
ple, Pjv calculated from the polarizability of the slab is 
P/v = 2a{N)Eioc/NV,B- CToo can be calculated from the 
bulk dielectric constant. The electrostatic of a macro- 
scopic slab in an external field Eioc gives: 



1 es-l 

Coo — 

47r es 
and from Eq. El we get: 



Eio 



a{N) = N 



Stt 



Poo'S' 

P, 



loc 



(19) 



(20) 



Therefore, the polarizability of an N layers slab is lin- 
ear in the number of layers, the slope of the straight line 
depends only on the bulk dielectric constant, and the in- 
tercept at the origin measures the difference between the 
dielectric properties of the bulk and of the slab surfaces. 
The fit gives es = 2.91 close to the value of the dielectric 



constant calculated in bulk benzene 



2.87. The term 



PooS/Eioc has the dimensions of a polarizability. Defin- 
ing as = ' "^here the factor 2 accounts for two slab 
surfaces, we olatain as = 4.5 a.u.. Therefore the surface 
layers are slightly more polarizable than the bulk, but in 
benzene the effect is indeed very small. Of course, this 
conclusion is valid only for the electronic contribution to 
the dielectric constant, and is obtained in a model sys- 
tem where the truncated bulk geometry is used for the 
surfaces. Finally, we note that, in benzene, the bulk unit 
cell contains two layers, but the polarizability of the slabs 
does not show any even-odd effect in Fig.Qfor symmetry 
reasons. 
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